library(hydroGOF)
library(clusterSim)

data1=read.csv("C:/data/mod_daily_merge.csv")
data1=subset(data1,data1$station!=3)

pm1=data1$pm1
pm25=data1$pm25
pm10=data1$pm10

aod=data.Normalization(data1$aod,type="n5",normalization="column")
temp=data.Normalization(data1$temp,type="n5",normalization="column")
avg_temp=data.Normalization(data1$avg_temp,type="n5",normalization="column")
avg_rh=data.Normalization(data1$avg_rh,type="n5",normalization="column")
avg_preci_24=data.Normalization(data1$avg_preci_24,type="n5",normalization="column")

data1=data.frame(pm1,pm25,pm10,aod,temp,avg_temp,avg_rh,avg_preci_24)
model1=lm(data1$pm25~data1$aod+data1$temp+data1$avg_temp+data1$avg_rh+data1$avg_preci_24)

pm25model=model1$fitted.values
model1_re=sum(abs(data1$pm25-pm25model)/data1$pm25)/nrow(data1)*100
model1_rmse=rmse(data1$pm25,pm25model)
model1_r=cor(data1$pm25,pm25model)

print("Model 1:")
print(paste("Number of sample: ",nrow(data1)))
print(paste("R2:",round(model1_r*model1_r,3)))
print(paste("RMSE:",round(model1_rmse,3)))
print(paste("RE:",round(model1_re,3)))


data1=subset(data1,abs(data1$pm25-pm25model)<3*model1_rmse)
model1=lm(data1$pm25~data1$aod+data1$temp+data1$avg_temp+data1$avg_rh+data1$avg_preci_24)
#model1=lm(data1$pm25~data1$aod+data1$temp+as.numeric(data1$area_id))
pm25model=model1$fitted.values
model1_re=sum(abs(data1$pm25-pm25model)/data1$pm25)/nrow(data1)*100
model1_rmse=rmse(data1$pm25,pm25model)
model1_r=cor(data1$pm25,pm25model)


print("Model 1:")
print(paste("Number of sample: ",nrow(data1)))
print(paste("R2:",round(model1_r*model1_r,3)))
print(paste("RMSE:",round(model1_rmse,3)))
print(paste("RE:",round(model1_re,3)))


